Auxiliary field Monte- Carlo for charged particles 
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This article describes Monte-Carlo algorithms for charged systems using constrained updates for the elec- 
tric field. The method is generalized to treat inhomogeneous dielectric media, electrolytes via the Poisson- 
Boltzmann equation and considers the problem of charge and current interpolation for off lattice models. We 
emphasize the differences between this algorithm and methods based on the electrostatic potential, calculated 
from the Poisson equation. 



I. INTRODUCTION 

Fast methods for calculating Coulomb interactions are of 
the greatest importance in the simulation of charged con- 
densed matter systems. Several algorithms are available in- 
cluding multipole' and Fourier based^ methods; they permit 
efficient molecular dynamics simulations but are ill-adapted 
to Monte-Carlo simulation. These methods have been widely 
used and implemented; recent research has led to multi-time 
step algorithms and aims at generalizing the methods to multi- 
processor environments. Despite this effort the Coulomb loop 
is still the slowest part of many simulation codes^. 

An alternative approach, the subject of this paper, is to use 
an auxiliary field in order to calculate the Coulomb interaction 
in a manner familiar from electromagnetism. Many books'*-^ 
consider the following scalar functional: 



d^r 



(1) 



where p is the charge density imposed by the experimental 
geometry and is a freely variable field. One might hope 
that one could use this functional to construct alternative, and 
more efficient algorithms due to the local nature of the energy 
function. One is instructed to minimize Eq. Q by taking the 
functional derivative with respect to cj); this leads to Poisson's 
equation familiar from electrostatics: 



(2) 



Thus we have a minimum principle; we can identify the func- 
tion (f)p which minimizes Eq. Q with the electrostatic poten- 
tial. However when we substitute the solution of Eq. (|2j in 
Eq. ([l) (using appropriate boundary conditions) we find 
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grad I 
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This is the negative of the true electrostatic energy. Recently^ 
a series of such potential functionals have been introduced in 
order to treat inhomogeneous dielectric systems, such as in- 
terfaces and pores, however all the functionals have a similar 
sign defect. 

This sign change is a major problem for numerical work. 
It implies that the same energy can not be used to evolve the 
field and the particle degrees of freedom in parallel. Clearly if 
one tries to do Monte-Carlo with cj) one generates the partition 



function 
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as can be seen by substituting (f) — (f>p + cj). This is the wrong 
statistical weight for particles interacting via Coulomb's law. 
While evaluation of the energy seems to have been reduced to 
a local calculation in practice a global optimization is needed 
every time step to minimize Eq. Q and then flip the sign of 
the result. For Monte-Carlo this leads to such a level of inef- 
ficiency that the method is useless. 

This sign problem is most surprising: We are aware that 
Maxwell's equations allow the co-evolution of field and parti- 
cles with the correct sign of the electromagnetic energy. What 
has gone wrong in the formulation of electrostatics in terms 
of an auxiliary field? The solution to the paradox is found by 
recognizing that the correct sign of interaction is only possi- 
ble with a vector field, a point made forcefully in Ref.^ in a 
general discussion of the forces of nature. 



II. VECTOR FUNCTIONALS 

Recently^ an efficient Monte-Carlo algorithm for charged 
systems was formulated using electric field lines (rather than 
the electrostatic potential) as the true dynamic degree of free- 
dom. The algorithm consists of a local energy functional 
together with local update rules; no global optimization is 
needed. Coulomb interactions result from the energy 



E 



(5) 



where the electric field, E is constrained by Gauss' law, 
div E — p/eo — . Introduction of a Lagrange multiplier 
for the constraint implies that the minimum energy is indeed 

l^p{4'p)^ J y(grad0j,)2d3r, (6) 

with (f)p (now the Lagrange multiplier) again the solution to 
Poisson's equation. 

We now show that the use of the energy Eq. (|5} at finite 
temperature still leads to pure Coulomb interactions. Let us 
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define a partial partition function by integrating over the elec- 
tric field but not the particle positions, r^: 

Z(r,) = /" PE <5( div E - p(r,)/eo) e''^" . (7) 

r 

Let us evaluate the integral by changing variables, Et^ = E + 
grad (pp. We find that 

Z(r,) = / VEtr n '5(divEt,)e-'^^.''(^"-s-d0p)^d^r 

r 

/ VEtr n^(divEt.)e-^*-^^-'^'- 

r 

= 2co«iom&(ri) X const 

The cross term in the energy, / Et^ • grad (j)p d^r, is identi- 
cally zero as can be seen by integration by parts. We discover 
that integration over the tranverse modes multiplies the statis- 
tical weight of each configuration by a constant. The longitu- 
dinal component of the electric field, — grad (f>p, produces an 
effective Coulomb interaction between the particles. 

We shall now give a full description of the algorithm show- 
ing how configurations can be dynamically generated accord- 
ing to the constrained statistical weight in Eq. (Q. This gives 
rise to a lattice gas in which charges are confined to a grid. 
We demonstrate numerically the efficiency of the algorithm 
for a system of lattice dipoles. A coarse graining of the elec- 
tric degrees of freedom leads to a new algorithm for treating 
the Poisson-Boltzmann equation as well as a non-constrained 
algorithm for the Yukawa potential. We then show how to 
treat off lattice models by interpolating charges to the grid. 
We show that the generalization of the algorithm to inhomo- 
geneous dielectric media automatically re-sums a large class 
of fluctuation based potentials such as the Keesom potential, 
capturing important characteristics of polar dielectric materi- 
als which are missed by the Poisson equation. 

III. PERIODIC BOUNDARY CONDITIONS 

1. Conventional approach 

Imposition of periodic boundary conditions in systems in- 
teracting with Coulomb interactions is subtle due to the con- 
ditional convergence of Coulomb sums'*: When extrapolating 
periodic images of the system to infinity there is a residual 
dependency on the dielectric constant, eg of the surrounding 
medium. 

The total potential, $ is the sum of (j)p, coming from the 
solution of the Poisson equation with periodic boundary con- 
ditions, and a dipolar contribution: $ = 0p + 0^- Here, 

0d = -ET, (9) 

where E is a vector proportional to the dipole moment per unit 
volume of the simulation cell, d 

1 



E 



eo(l + 2e,) 



(10) 



This constant field reminds us of the "depolarizing field" in el- 
ementary treatments of spherical dielectric mediai. Two com- 
mon choices of boundary condition are "tin-foil" with eg = oo 
and "vacuum", e, = 1. Clearly the tin-foil limit corresponds 
to keeping just the periodic potential, (f>p. 

The field E generates a purely additive contribution to the 
energy density: 
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This extra energy term can be awkward to handle in molecular 
dynamics; its definition leads to discontinuities in the energy 
function since = CiVi is calculated with a given, 

fixed Bravais lattice; V is the simulation volume, the charge 
of particle i at position in the simulation cell. 



2. Constrained algorithm 



The general solution to Gauss' law in periodic systems is 



E = - grad i 



curl Q + E 



(12) 



where Et^ — curl Q is an arbitrary transverse field. 
Our algorithm is closely related to the following Maxwell 
equation:"' 



dt 



curl H , 



(13) 



with J the electric current and H the magnetic strength. If we 
start a simulation in a state in which E = then an evolution 
of the field according to Eq. ( I13> (as well, it will be seen, as 
with our algorithm) generates fields with 



dt 



d^r J/eo -d/eo 



(14) 



so defining d. We have used the fact that the integral of the 
curl of a periodic function is zero. This is very similar to 
the Ewald convention of Eq. ilQi with Eg = 0. We note that 
this choice is rather unconventional, one normally considers 

1 < Eg < OO. 

If the charges in the simulation are bound (so as to remain 
in the basic unit cell) then d = d^;. If there are free charges 
in the simulation, ±e, then the same configuration can be pro- 
duced by different currents which wind about the periodic sys- 
tem: We only require that d — dg = eh/V, where b is a 
Bravais lattice vector We again find that there is an additive 
contribution to the energy density of the form 



Uo 



d^ 

2eo 



(15) 



We note the natural use of d rather than d^ is a great advan- 
tage for molecular dynamic simulation. It leads to energies 
which are continuous functions of the particle position. 

In our simulations of periodic systems we are also capa- 
ble of reproducing tin-foil boundary conditions. We do this 
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by adding a constant electric field Eg to the zero wavevector 
component of the electric field, so that E — — d/eo + Eg. 
This leads to a simple modification of uq, 



Uo 



2eo 



(16) 



Integration over Eg as well as the transverse field now 
leads to an effective statistical weight which varies as 
g-/3/2/(grad ^pf d^r^ particlcs interact through the pe- 
riodic solution to Poisson's equation, (j)p. 

Finally in this section we note that a consistent description 
of polarization effects in quantum mechanics (based on con- 
siderations of gauge invarience) require the use of d rather 
than d e for the polarization" . 



IV. A SCALAR FUNCTIONAL VIA GAUGE 
TRANSFORMATIONS 

In this section we continue the development of our formula- 
tion of the electrostatic potential in order to show that despite 
the failure of the variational principle Eq. Q it is possible to 
link the electrostatic energy Eq. (|5|i to a scalar potential. This 
section is independent of the practical implementation of the 
algorithm and may be omitted on first reading. 

Our separation of the electric field into longitudinal and 
transverse components is similar in many ways to the 
Coulomb gauge in conventional electromagnetism. It leads, 
however, to a description of the dynamics in terms of poten- 
tials which is non-local; as in conventional electromagnetism 
other choices of gauge are possible. 

Let us write the electric field in an arbitrary gauge in the 
form 



E 



■ grad ( 



From Gauss' law we deduce that p/eo 
now apply a gauge condition 



div A + a— = 
ot 



(17) 
div A. We 

(18) 



for some constant a. We discover that the propagation of the 
potential </) is now local: 



"1^ = V20 + p/eo 
ot 



(19) 



The dynamics of the electric field in the algorithm are related 
to the Langevin equation'" 

9E 

— = (eoV^E - grad p) /i - J/eq + C(0 (20) 

where ^ sets the time scale of the dynamics and C,{t) is a 
transverse noise. From these equations we see that the spe- 
cial choice e^a = ^ allows a particularly simple equation for 
the vector potential A. 



dA _ 1 
dt a 



V^A + J/eo + CW 



(21) 



In this diffusive gauge we calculate the electric energy Eq. Q 
in terms of the potentials 



^2 . I d^r (22) 



Cross terms in A • grad (p have been eliminated with help 
of the gauge condition and by using the equation of motion 
for the scalar potential. We recognize that this functional is 
closely related to the naive scalar functional Eq. Q; it is, how- 
ever, dynamically constrained by Eq. (I18> 



V. ALGORITHM 

The Monte-Carlo algorithm uses the Metropolis criterion 
together with the constrained energy function Eq. (|5jl. In order 
to generate configurations according to the statistical weight 
of Eq. we need to 

• move particles without violating the constraint of 
Gauss' law to preserve the delta function on div E — 

• integrate over the transverse degrees of freedom Et^, of 
the electric field 

• integrate over Eg. 

Our third Monte-Carlo step eliminates the slaving of E to 
the current. If we perform simulations in which we do not 
integrate over the variable Eg we find a summation of the 
Coulomb energy which is different from that generally used 
in the Ewald method, but identical to that implied when inte- 
grating the Maxwell equations. 



3. Implementation 

The system is discretized by placing charged particles on 
the N = vertices of a periodic cubic lattice, {i}. The 
components of the electric field Eij are associated with the 
3N links {i,j} of the lattice. There are 3N plaquettes on the 
lattice each defined by four links. We use the notation Ei_2 to 
denote a local contribution to electric flux leaving 1 towards 
2. If the link from i to j is in the positive x direction we 
consider that this is the local value of the x component of the 
field E. The 3N variables Eij are thus grouped into TV three 
dimensional vectors. 

It is convenient to consider Gauss' law in the equivalent 
integral form 



(23) 



where is the total charge at the site i, enclosed by the surface 
integral, which we discretize as 



(24) 
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A is the lattice spacing. The sum in Eq. i24\ corresponds to 
the total electric flux leaving the site i through the surrounding 
plaquettes. Since we are considering Eij as a directed flux we 
have Eij = ^Ej_i. 

In order to integrate over E we add a global background 
vector Eg to Eij . The background field does not contribute to 
the divergence of the electric field, it does, however contribute 
to the energy 

U=^Y.(E.-^+Kf^ (25) 

links 

with Eg the component of the background field parallel to the 
link i, j. For convenience we also store the sum of the vari- 
ables Ei j, we denote this single 3-vector 'Eitotai- For instance 
the X component, Ef^^^i — Ei j, where the sum is over 
all X directed Hnks. We note that E — Etotai/N + 'Eg. 

With this discretization of the electric field we check the 
derivation of the Coulomb interaction between charged parti- 
cles. Particles now interact via a lattice Green function, rather 
than with the continuum Green function of Gc = 1 /r: The 
interaction is identical to that found from solving the Poisson 
equation with a simple difference scheme on the same lattice. 
Explicitly we find that 

' J (27r)3 6 — 2cosg^ — 2cosgy — 2cos(72 ^^^^ 

The lattice also regularizes the self energy of the particles 
which scales as ef / eA. This self energy is analogous to the 
Born self energy important in solvation theory. 

4. Particle Motion 

The simulation starts with Gauss' law satisfied as an initial 
condition. We displace a charge, e situated on the lattice site, 
1, to the neighboring site, 2; for figures and more detail we 
refer the reader to our previous work^. The constraint is again 
satisfied after motion of the particle if the field associated with 
the connecting link is updated according to the rule i?i 2 
E1.2 — e/eoA^. The trial is accepted or rejected according to 
the Metropolis algorithm. If the move is accepted we must 
remember to modify Etotai- It is this update that leads to the 
slaving of E to the current, Eq. il3\ . since it corresponds to 
eoE = -J 



5. Plaquettes 

We integrate over all transverse modes in the partition 
Eq. by modifying together the four link field values of one 
of the 3A^ plaquettes. In such an update the sum of the en- 
tering and leaving fluxes at each site does not change. This 
update modifies the electric field by a pure circulation; Etotai 
does not change. We perform this update by choosing one of 
the plaquettes randomly. The amplitude of the trial up- 
date is uniformly distributed between —Of) and 9q where 9q is 
chosen to have an acceptance rate of between 40% and 60%. 



Let us note that the distinction between link and plaquette 
degrees of freedom is very similar to that found in the Yee 
algorithm for integrating Maxwell's equations'^. 

6. Background field 

The Metropolis algorithm for the variable Eg is assured by 
using the energy Eq. J25> . However, since the individual link 
fields, Eij, are unchanged when we modify Eg we calculate 
the change in energy using the simplified expression 

^9 = ^ (NEl + 2Eg . Etotai ) . (27) 

By following the evolution of Etotai we avoid having to cal- 
culate the sum of the electric fields; modification of the global 
field is possible in a time which is 0(1). 

VI. EFFICIENCY 

In— we gave an demonstration of the efficiency of the al- 
gorithm for a charged lattice gas; in this section we hope to 
convince the reader that the algorithm is equally efficient in 
simulating dipolar systems in the absence of Debye screen- 
ing. 

To model a dipolar system we take a lattice gas of charged 
(and self-avoiding) particles as described above and pair each 
positive charge with a negative charge: They are linked by a 
harmonic, zero length, spring so that the elastic energy be- 
tween them is given by Ugpring = 7(rj ~ rj)^/2, is the 
position of the particle i on the lattice, 7 the spring stiffness. 
During the simulation we record the Fourier components of 
several important quantities: the charge density, the particle 
density and the transverse electric field (that part of the field 
for which the discrete divergence Eq. ( I24> is zero). From 
these recordings we deduce the two-time correlation func- 
tions. We find that the temporal correlations of the Fourier 
components are well described by single exponentials, as de- 
scribed previousljii^. 

We now plot the decay rate of each mode measured in 
"particle-sweeps" as a function of wavevector. Fig. Q. One 
particle-sweep is defined as the time needed to attempt (on 
average) one update for each particle in the simulation, in- 
dependent of the number of Monte-Carlo trials on other de- 
grees of freedom; we update the plaquettes with various rates 
in order to study the effect on the dynamics. Measuring time 
in this way allows us to study the intrinsic dynamics of the 
dimers and their coupling to the transverse field. Note that the 
computational effort needed to move a particle or to update a 
plaquette is almost the same; both are dominated by the time 
needed to generate random numbers and calculate the expo- 
nential function needed for the Metropolis algorithm. In our 
simulations we used a lattice of = 20'' sites and 24,000 
plaquettes with 1200 dimers (2400 particles). 

In looking at Fig. Q one should first note that there are 
two diffusive modes, with the decay rate proportional to q^, 
where q is the wavevector of the mode. They correspond to 
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density fluctuations and to electric field diffusion, described 
by Eq. j20t . The third mode, describing charge fluctuations, 
has a finite frequency at long wavelengths; it couples to the 
internal modes of the dimers which relax rapidly due to the 
spring. 
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FIG. 1: Decay rate, in particle-sweeps, of correlation functions as a 
function of wavevector. □: density fluctuations. A: charge fluctua- 
tions, O- transverse electric field. 60,000 particle-sweeps, 8 record- 
ings per sweep. A = 1, temperature T = .5, |e| = 1, 7 = 1/2, 
€= 1. The single mode at q = describes the dynamics of E. 

For illustrative purposes we have tuned the relative proba- 
bilities of particle motion and plaquette updates in Fig. Q so 
that the diffusion rates of the density and the transverse elec- 
tric field are very similar. To do this we used a relative prob- 
ability (per degree of freedom) of particle motion to plaquette 
updates of 1 : 1/3; Since there are considerably more plaque- 
ttes than particles about 75% of the time in the simulation is 
devoted to updating the transverse field. For simplicity we do 
not update Ej,; we are simulating with "Maxwell" boundary 
conditions. 

In Fig. (13 we simulate a similar system with the same ratio 
of particle motion to plaquette updates. The only modification 
is that the charges on the dimers are now put to zero; there is 
no charge in the simulation. We note two points in comparison 
with Fig. Q: The diffusion rate of the dimers has increased 
by some 60%, while the diffusion rate of the electric field has 
slightly fallen. In this plot we replace the charge-charge cor- 
relation function by an analogous partial structure factor, in 
which the scattering amplitude of one particle in a dimer is 
+1, while the second particle has a scattering amplitude of 
— 1. This structure factor is again sensitive to the internal dy- 
namics of the dimer. 

In Fig. Q we again simulate charged dimers, we reduce 
however the number of plaquette updates in the simulation so 
that they now occur in the relation 1 : 1/30. In this simulation 
75% of the time is spent updating the position of the particles, 
only 25% of CPU time is devoted to the transverse field. As 
might be expected the propagation of the transverse modes 
has slowed down. However the equilibration of the charge 
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FIG. 2: As Fig. Q but with charge-free dimers. Charge correla- 
tion function replaced by partial structure factor sensitive to internal 
dimer motion. The dimers are diffusing faster. 



and density modes is not greatly modified in comparison with 
Fig. Q. 
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FIG. 3: As Fig. Q but with fewer updates to transverse modes. 
Dimers carrying charges. The relaxation dynamics of density fluctu- 
ations and the charge fluctuations are only slightly modified despite 
the order of magnitude decrease in the work done on the plaquettes. 

In Fig. (0} we have simulated charge-free dimers with pla- 
quettes updated in the proportion of 1 : 1/30. We note that 
compared with Fig. Q (where the dimers carry charges) the 
transverse modes have substantially slowed. Motion of the 
particles is playing an important role in agitating the field vari- 
ables in Fig. (|3j; this accelerates the propagation of the elec- 
tric field. 

From our results on the simulation on this simple dimer 
system we find that quite remarkably one can "get away" with 
very little work on the plaquette degrees of freedom. How can 
this be so? Motion of a particle along the four links defin- 
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ing a plaquette is equivalent to a single plaquette update; the 
motion of the particles is, on its own, quite good at integrat- 
ing over the transverse degrees of freedom. This integration 
is clearly imperfect since the only values of the plaquette cir- 
culation that are generated are integral multiples of e/eA. We 
see, however, that the plaquette updates only have to "fill in" 
between the integral multiples of e/eA rather than perform the 
full integration of these modes. 

It is also rather natural that longitudinal and transverse 
modes (with scalar and vector symmetries) are weakly cou- 
pled in a homogeneous system at long wavelengths; the low- 
est order couplings occurring in the free energy are presum- 
ably at least cubic; if the number density fluctuation is n and 
the transverse field fluctuation etr one could expect that the 
lowest order coupling in a Landau picture is e^^n. 
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FIG. 4: As Fig. Q but with fewer updates to transverse modes 
and charge-free dimers. The dynamics of the transverse modes have 
slowed in comparison with Fig. J3} in the absence of coupling to the 
charges. 



Examining the figures we see that the total overhead in sim- 
ulating the charged dipole, compared to a charge-free dimer is 
about a factor 2 with appropriately chosen simulation param- 
eters. 

The simulation of Fig. Q corresponds to 30 minutes of 
simulation on a pentium 4 2.5Ghz. Total clock time was 
(much) longer due to the overhead of the Fourier analysis of 
the mode structure. 

We note that the temperatures used in these simulations is 
very high compared with those of interest in condensed mat- 
ter physics. The Bjerrum length is comparable to the particle 
size in these simulations at a temperature of about T = 1 /An, 
the present simulations are however at T — 1/2. We find 
that at low temperatures the mobility of the particles drops 
very strongly, since there is a finite barrier for particle motion 
due to the discrete nature of the charges and the lattice. We 
have recently implemented an off lattice version of the algo- 
rithm and find that we can easily overcome this artefact of the 
discretization'-^. 



VII. POISSON-BOLTZMANN SIMULATIONS 

Classic methods for treating the Poisson-Boltzmann equa- 
tion require a calculation of the long ranged Coulomb interac- 
tion, either in real space or in Fourier space''*. In this section 
we show how constrained Monte-Carlo gives rise to a local 
treatment of the Poisson-Boltzmann equation. The method, 
however, has some original features since rather than consid- 
ering the minimum of the Poisson-Boltzmann functional we 
integrate over the charge degrees of freedom, creating a fluc- 
tuation enhanced Poisson-Boltzmann algorithmic^ 



A. Formulation of the free energy 

Consider a system of n particles in a fixed volume V with 
coordinates q. The partition function is 



dq, 



(28) 



where U is the configurational contribution to the energy and 
the integral is over all particle positions. Let us coarse grain 
the evaluation of the partition function by subdividing the vol- 
ume in to a large number of elementary cells of volume A'^, 
replacing the configuration integrals by sums. In each cell let 
there be rn particles. Then 



-f3U 



(29) 



Here the combinatorial factor counts the number of ways of 
distributing particles between boxes and the energy has been 
re-expressed in a coarse grained manner, U = U{ni), in terms 
of the cell occupation numbers. The factor A"^"' comes from 
the rii fold integral over the cell volume. Re-expressing the 
rii in terms of local concentrations of particles, pi — n^/A^, 
and exponentiating the combinatorial factor (using log(n!) « 
nlogn — n ) we find the effective energy function for the 
coarse grained description 

H= fcBrA3^(pan(p,A3)-p,)+z:/({pJ), (30) 

i 

where the sum over i is now a sum over both boxes and chem- 
ical species. This expression is very similar to that used in 
many mean field studies, but the derivation shows that this 
effective energy function is valid for fluctuating densities. It 
is not a functional simply valid as a minimum principle as is 
emphasized in density functional approaches to the Poisson- 
Boltzmann equation. In our simulations the variables pi are 
fluctuating variables. 

In order to describe a charged coarse grained fluid we 
must choose A: The ions are separated on average by a 
distance ^ — p^^l'^ while the Bjerrum length is given by 
l\, — e^/AiT ksT e. We also introduce a further length scale: 
the Debye length, the length scale at which the energy stored 
in the electric field, due to concentration fluctuations, is equal 
to ksT ; we find that £^ ^ £,^/^b- In order to use the coarse 
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grained entropy we require that A ^ ^ so that ^ 1, how- 
ever if we wish to describe screening in detail we also require 
K ^ td implying a relatively narrow window of possible val- 
ues for the coarse graining parameter 
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FIG. 5: A positively charged plane is simulated in in the presence of 
counter-ions using the Poisson-Boltzmann formulation ea J31> . We 
plot the concentration of the of two, synmietric ion concentrations as 
a function of the distance from the plane. System of size 50 x 50 x 
50. Concentration is for a single equilibrated system. A negative 
charges, positive charges. 



The full system we simulate is 



(31) 



i 

A3 ^ 2 

links 



where the electric field is constrained by Gauss' law, with both 
external charges and pi acting as sources. 

Other contributions to the free energy varying in p3/2 ^j-g 
often introduced as local field corrections in density functional 
theories. However, such energies are derived on a length scale 
larger than the Debye length and are due to correlations of the 
electric field over id- In the effective free energy measured on 
a scale A < id such terms should not be added by hand. They 
will be generated automatically by the algorithm via fluctua- 
tions in the charges. 

The derivation of the effective Hamiltonian on the scale A 
is rather close in spirit to the renormalization group. It would 
be interesting to perform the numerical renormalization group 
on a charged fluid in order explicitly measure the effective 
Hamiltonian as a function of the scale for £^ < A < id- The 
derivation breaks down in some particularly interesting cases, 
including strongly charged surfaces. In this case a new, inde- 
pendent, length, the Gouy-Chapman length becomes smaller 
than the coarse graining cell. 

As a test of the algorithm we simulated. Fig. (|5}, a uniform 
positively charged plane in a background of positive and neg- 
ative ions. We see the depletion of the positive ions and the 



attraction of negative ions near the surface. In the simulation 
the plane was stationary however the method allows parallel 
dynamics of both the ionic and the macro-ion degrees of free- 
dom without having to perform global updates or optimiza- 
tions. Such simulations are useful in the simulation of flexible 
charged objects such as polyelectrolytes within the Poisson- 
Boltzmann approximation; in such flexible objects we can 
individually update the position of individual charges in the 
macro-ion as described in the application to electrolytes. The 
method is less likely to be useful in the simulation of large 
rigid objects such as a colloidal sphere carrying large charges. 
In such cases the simultaneous modification of E on many 
links leads to a low Monte-Carlo acceptance rate, or when 
working off-lattice small step sizes. Our method is also in- 
compatible with smart global algorithms, such as the pivot 
algorithm in polymer simulation, for similar reasons. 

In the functional Eq. ( 13 U we integrate over the continuous 
variables pi. However, numerically it is just as easy to keep 
the discrete charges, n^. In such a way one still describes a 
coarse grained system so that (ni) > 1 but a higher degree of 
detail is retained. A improved account for the entropy is also 
possible if a better approximation for log (n! ) is used. It would 
be particularly interesting to understand if attractive, correla- 
tion effects can be reproduced within such a coarse grained, 
discrete model. 

Finally we wish to make some remarks about the ergodicity 
of the Poisson-Boltzmann algorithm: We have already seen in 
our dipolar model that the motion of the charges excites trans- 
verse degrees of freedom of the field. It was only the discreet 
nature of the charges in the dipolar system which prevented 
the transverse field from coming to equilibrium even in the ab- 
sence of independent field updates. In the Poisson-Boltzmann 
algorithm the amount of charge transfered at each step is ran- 
domly chosen: the algorithm is ergodic even without the field 
updates which can be dropped from the algorithm. 



B. Screened Coulomb Interactions 

In even coarser-grained treatments of charged systems the 
counter ion degrees of freedom are eliminated entirely by us- 
ing an effective Yukawa interaction 6162 exp(— Kr)/r. One 
implements this potential numerically by introducing a non- 
constrained energy function for the electric field: 



Uy 



fo 
2 



^(divE-p/eo)'] 



(32) 



The analytic treatment of this energy is simpler than the con- 
strained functions that we have considered until now; the 
functional is Gaussian in the field E. Clearly the limit k ^ 
leads to a constrained electric field. Taking the functional 
derivative of Eq. ( I32> with respect to the electric field we find 



^E — grad div E = — grad p/eo 



(33) 



for the stationary points of the energy function. This equation 
is best treated by separating the transverse and longitudinal 
components of E. The longitudinal component of the electric 
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field gives rise to the Yukawa potential at the minimum of 
Eq. i32\ . The transverse components of E do not couple to 
the density. 

When K is small the screening length becomes large and 
the field becomes more and more constrained by Gauss' law. 
Freely chosen Monte-Carlo moves with Eq. (13 2> then become 
inefficient since small amplitude trials are needed. It becomes 
more efficient to separate the Monte-Carlo moves into two 
classes, those which conserve div E — p/ep and those which 
modify it. 



Vm. OFF LATTICE INTERPOLATION AND RELATION 
TO PLASMA SIMULATIONS 

We now formulate the algorithm when the particles are 
moving in the continuum but with charges interpolated to the 
cubic grid. We shall see that the simplest generalization of 
our method is closely related to numerical methods used in 
plasma simulation. 

Consider a single particle within a cell of the grid. We here 
take A = 1 and consider a cube such that < {x, y, z} < 1. 
The simplest interpolation of the charges to the lattice is a 
linear interpolation so that the charge at the site (0, 0, 0) is 
equal to e(l — x)(l — y)(l — z). Interpolating to the eight 
sites of the enclosing cube preserves both the total charge and 
the dipole moment. 

When the particle moves the rate of change of the density 
at the origin is 

= -e[(l-x)(l-y)i + il~x){l-z)y 

+ (l-y)(l-z)i].(34) 

This variation has a simple, direct interpretation as currents in 
the links: 

J,=eil-x){l-y)z (35) 

is the current in the z direction on the link between (0, 0, 0) 
and (0,0,1). 

Eq. ( I34> is a perfect differential so that when the particle is 
transported around a loop 



F{t)dt = , 



(36) 



implying that the interpolated charge returns to its initial 
value; this is a clearly needed consistency requirement for the 
current. 

Consider now Monte-Carlo: We need to generalize the 
expression for Jz to finite displacements, while preserving 
Eq. i36\ . Take a path that starts at {xi,yi,Zi) and finishes 
at (xi+i, Zi+i) after a time t = 1. We integrate the ex- 
pression, Eq. (I35> for the link current to find 



Ae,/e = f 
Jo 



{1 ~ Xi - Vj:t){l - J/j - Vyt)Vz dt 



(1 - Xi){l ~ yi)vz 

{VxiX - yi) + Vyjl - Xi))Vz 

2 

VxVyVz 



where Vx = Xi+i — Xi. By construction a series of such steps 
satisfies an expression equivalent to Eq. i36\ . The transfered 
charge Ae^ leads to a corresponding modification of the elec- 
tric field on the link E E ~ Acz/eo. 

While being a low order interpolation scheme Villasenor 
and Bunemani^ showed that such a Gauss' law conserving 
scheme allows long time, stable integration of plasma equa- 
tions. Many other (non-conserving) methods need a peri- 
odic correction step in which Poisson's equation is solved to 
re-initialize the longitudinal electric field, leading to compli- 
cated and slow code. We again see the central role played 
by Gauss' law in the efficient simulation of charged systems; 
codes based on Buneman's TRISTAN would seem to be the 
closest analogy in the literature to our algorithm. 

For the simplest, low precision modeling of polymers and 
polyelectrolytes we expect that this present scheme will al- 
ready be useful. The scheme needs to be improved by interpo- 
lation over several cells in order to be competitive in accuracy 
with those that are currently used in the soft condensed matter 
simulation community^. We present one possible algorithm in 
an accompanying paper'^ on charge and current interpolation. 



IX. INHOMOGENEOUS MEDIA 

The forces on charges in the presence of inhomogeneous 
dielectric media are complicated; a charge is attracted towards 
regions of high e. Such forces are well known to be important 
in the structuring of charge distribution around proteinsii and 
at surfaces. 

The electric energy of an arbitrary inhomogeneous medium 
is giveni^by 



D2 



2e(r)' 



(38) 



where the electric displacement, D = e(i")E, obeys the con- 
straint div T) — p . The generalized Poisson equation found 
from this constrained variational problem is 



div (e(r) grad (fy) — —p, 
so that the electrostatic energy is given by 

/cfr) 
(i^r— (grad 0^ 



(39) 



(40) 



The solution to the constraint equation in non-periodic media 
is 



(41) 



(37) 



D = — e(r) grad (f> + curl Q . 



As for the the case of homogeneous media the energy Eq. ( I38t 
can be written as a sum of two terms corresponding to elec- 
trostatic and fluctuation contributions. 

The statistical mechanical treatment of inhomogeneous di- 
electric media is however more complicated than the case of 
homogeneous systems: We have seen that when e is a con- 
stant the contribution of the transverse fluctuations to the free 
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energy is independent of the positions of the particles; the al- 
gorithm gives rise to pure Coulomb interactions between par- 
ticles. Here we shall show that in addition to the electrostatic 
interactions, Eq. ( I40> . we also generate extra dipole-dipole in- 
teractions. 

We shall show that the fluctuation potentials are expected 
on general grounds in inhomogeneous media, they correspond 
to the Keesom potential''* varying as between fluctuating 
permanent dipoles; These interactions should be distinguished 
from the van der Waals interaction (also in 1/ r^) which has its 
origins in quantum mechanics^ 



A. Dipolar forces 

In this section we consider the fluctuation potentials in the 
absence of charges, then: 



n^(divD(r)). (42) 



This energy is a function of the position of the particles due to 
the presence of e(r) in the denominator of the energy. 

We will proceed by firstly re-writing the constrained field 
D in terms of a projection operator This is easily done in 
Fourier space where one can trivially separate longitudinal 
and transverse components of a vector field. We then per- 
form an inverse Fourier transform to find the analytic form of 
the projection operator in real space. Rather surprisingly we 
find that this projection operator is closely related to dipolar 
interactions which vary in An expansion to second or- 

der in the fluctuations in the dielectric constant then gives a 
contribution to the free energy which naturally decreases as 

One calculates a transverse field from a general uncon- 
strained field D by subtracting off the longitudinal compo- 
nent: 



Dt(q) =D(q) 



q(q.D(q)) 



(43) 



Since we extract the transverse component by multiplication 
in Fourier space the same operation can be written as a con- 
volution in real space: 



Dt(r) = J d^n r(r,ri)D(ri; 



(44) 



with a real space projection operator T. This inverse Fourier 
transform is a little tricky to perform: The projection oper- 
ator does not decay "nicely" for large q. We proceed more 
carefully using various operator identities: 

We consider the function v(r) • Dt(r) for an arbitrary vec- 
tor v(r). We note that iq-, iq and l/q^ are respectively the 
Fourier transforms of div , grad and 1 /Anr. This results in 



v(l)-Dt(l) 



v(l)U(l,2) 



1 



47rri; 



■grad divD(2)j d2 
(45) 



We have used a common shorthand notation replacing by 
the index i. We now integrate by parts twice (using the fact 
that the operators div and — grad are mutually adjoint) to 
transfer the derivatives from D to v/r while discarding sur- 
face terms: 

v(l)-Dt(l) = Jb{2) (^S{1,2)+ grad div ^-i— ^ v(l)d2 

(46) 

The operators " grad " and " div " act on the "2" degrees of 
freedom so that v(l) is to be considered as a constant. One 
simplifies using the identity 

grad (v • grad (l/r)) = (v • grad )(grad (l/r)) (47) 

3(v ■ f )f — V 



for r 7^ 



leading to 



(48) 



I is the unit matrix. The coefficient of the delta function is a 
little subtle but is found by noting that the trace of an opera- 
tor is independent of the basis: A projection operator which 
selects the two transverse degrees of freedom has trace two. 

The transverse projector is closely related to dipolar inter- 
action between particles. For instance the electrostatic inter- 
action between two dipoles pi and p2 when ri ^ r2 is just^« 



-piTi,2P2/eo 



(49) 



as is discussed in many standard texts'*. T has the useful prop- 
erty 



r(0,l) •T(l,2) d^l = T{0,2) 



(50) 



since it is a projection operator and thus idempotent. Here "•" 
signifies matrix multiplication of the operator T. 

We now proceed by perturbation theory in fluctuations in 
the dielectric constant. Let us consider two small inhomo- 
geneities in the dielectric properties separated by a large dis- 
tance; the background dielectric constant is eq- One expands 
the partition function to second order in the perturbation. Each 
fluctuation in the dielectric constant has associated with it an 
amplitude (l/e(r) — 1/eo) ~ —6e/eQ. 

The second order perturbation in the partition function due 
to interaction between two sites "a" and is found to be 



A2(a,6) 



44 



6e{a)6e{b) ^^{1234} 

X D(l)r(l, a) •r(a, 2)D(2) 

X D(3)T(3,6)-T(6,4)D(4) (51) 



which needs to be averaged with the statistical 
weight from Eq. (I42t . The interesting pairings are 
(D(1)D(4))(D(2)D(3)) and (D(1)D(3))(D(2)D(4)). 
To perform the integration over transverse modes we 
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introduce the modified energy functional 



B. Simulations 



U 



D2 

2^' 



D(q) = (l + A|q)(q|)D(q), 



(52) 
(53) 



In the limit of A large and positive all longitudinal fluctuations 
are suppressed in the measure, however we can now perform 
usual, unconstrained perturbation theory with this Gaussian 
energy. In the limit of large A the appropriately normalized 
correlators become 



(54) 



The calculation is much simplified by making use of the 
identity Eq. ( 15 Ot . The free energy cost of dielectric inhomo- 
geneity is given by 



F = - keT InZ = Fo + V{a, h) 



(55) 



where Fq is the free energy of a uniform dielectric medium. 
The potential between particles is found to be 



V{a,h){r) = - 



3 ksTVaSCa VbScb 



TTT^{a,b) (56) 



(57) 



a.b 



with Vi the volume of the dielectric inhomogeneity. Substi- 
tuting that the relative dielectric constant for a weakly dipolar 
material (with ei/eo — 1 small) is given by 



3eo ksT 



(58) 



with Hi the density of dipoles of strength pj we recognize 
the Keesom potential''' between pairs of thermally agitated 
permanent dipoles 



Vk 



9 9 
PaPfc 



3 ksT {Aneor 



3\2 



(59) 



Note the distinction between free energy and internal energy 
for the Keesom potential. The internal energy is given by 



(60) 



When the dipole moments are independent of the temperature 
UK — 2Vr-. If Ei is independent of the temperature uk — 0; 
the potential is purely entropic. 

The presence of these dipole terms seems inevitable in the 
Monte-Carlo algorithm: Removing them requires the evalua- 
tion of a complicated multibody determinant. The algorithm 
is closer in spirit to simulations which replace the dielectric 
medium by an ensemble of fluctuating Langevin dipoles^ 
than the classic, zero temperature, solution of the Poisson 
equation. 



Let us discretize the energy Eq. ( I38> by placing the particles 
on the vertices of a cubic network. To each lattice point we 
also associate a dielectric constant e,;. The energy associated 
with a link "i, j" is then defined to be 



1 



4 \e. 



1 



(61) 



The constraint J^j Dij = is imposed as an initial con- 
dition on D rather than E. The electric displacement D is 
updated in a similar manner to that described above for the 
electric field. The algorithm is simple compared with direct 
solution of the Poisson Eq. ( I39> ; in particular we note that 
Fourier methods are do not diagonalize the Poisson equation 
when e is a function of the position. 
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FIG. 6: Plot of the density-density structure factor as a function of 
wave-vector. 1000 neutral particles simulated on a 15 x 15 x 15 
grid, o: ep = eo- *: = 5eo- °' ~ 0.2eo. For small q the 
structure factor increases, reflecting attractive interactions between 
the particles when 7^ eo- 

We implemented the algorithm for a neutral lattice gas of 
dielectric constant ep embedded in a background eo in the ab- 
sence of free charges. The density-density correlation func- 
tions S{q) were measured to characterize the effect of interac- 
tions. Fig. (|6j. As expected when gp = eo no structure appears 
in the simulation. For both Ep < eq and Ep > eq the structure 
factor increases for small q. This can be interpreted as being 
due to an attractive interaction between the particles leading 
to dumpiness in the structure. 



C. Low symmetry materials 

In low symmetry materials the dielectric constant is a ten- 
sor, with up to three distinct eigenvalues and eigenvectors. 
The algorithm generalizes to this case; in such materials the 
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constraint equation is 

div D = p , (62) 

where D = eE; e is a matrix to that D and E need no-longer 
be parallel. The electric energy 

^=\j De^^D d=^r. (63) 

denotes the matrix inverse. 
One possible application of this discretization is in liquid 
crystals where the order parameter manifests itself as a dielec- 
tric anisotropy. Models are easily constructed within a mean 
field description. This permits numerical study of transitions 
driven by electric fields. 

X. CONCLUSIONS 

This paper has presented local algorithms for the simulation 
of charged systems. Our methods are strongly inspired by the 
observation that Maxwell's equations allow a local evolution 
of field degrees of freedom leading to the Coulomb interac- 
tion. It introduced dynamics which conserve div E — p/e 



while starting the simulation with this constraint satisfied. The 
Monte-Carlo algorithm eliminates the magnetic degrees of 
freedom from Maxwell's equations leading to pure electro- 
static interactions. 

The method is rather general and can be used in many situa- 
tions: for instance application of the constrained formalism to 
path integral Monte-Carlo should be equally straight forward 
as the cases treated here. The main weaknesses that we see in 
the algorithm are 

• Incompatibility of constraint preservation with simula- 
tions in the grand-canonical ensemble, so useful in the 
study of phase behavior. 

• Low efficiency with large strongly charged objects un- 
dergoing coherent motion, either due to rigidity or the 
use of smart global moves. 

The code used in the simulations is available from the au- 
thor 

1 would Uke to thank Ralf Everaers for many discussions, 
essential in the formulation of these algorithms. 
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